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Abstract 


Introduction. Mathematical models and methods are widely used to study natural phenomena, replacing more 
expensive field experiments. However, one of the main challenges in modeling processes in complex systems is the lack 
of available input data and difficulty in selecting model parameters. The use of observational data assimilation methods 
is one of the ways to provide mathematical models with input data and parameter values. The aim of this study was to 
predict the development of complex natural systems under conditions of pollution using mathematical modeling 
techniques. To achieve this, several tasks were completed: a method for assimilating observational data was selected, a 
mathematical model for biological kinetics was updated, it was integrated with a hydrodynamic model, and a software 
package was developed. The significance of the work lies in the to the implementation of a model of the dynamics of 
phytoplankton populations (eutrophication) of the Azov Sea in the presence of pollutants, based on the use of 
variational methods for assimilating data obtained during expeditionary research. 

Materials and Methods. The spread of pollutants was modeled using a three-dimensional mathematical model based on a 
system of convection — diffusion — reaction equations. The vector of movement of the aquatic environment was the input 
data for the model. The components of the current velocity vector in the coastal system were calculated using a mathematical 
model of hydrodynamics, based on three equations of motion and the equation of continuity. The software package developed 
based on these models received full-scale data collected during expeditionary research as input, and allowed us to refine the 
model of pollution in the aquatic environment and biota using variational methods for data assimilation. 

Results. A short-term forecast for the spread of pollutants at the outlet of the Taganrog Bay was developed. The 
conducted computational experiment reflected the dynamics of pollutant spread from sources of contamination over a 
period of 3 to 12 days. 

Discussion and Conclusion. The variational methods of assimilating observational data discussed in this study allow 
for the refinement and supplementation of mathematical models of phytoplankton population dynamics and pollutant 
spread. The software based on these mathematical models enables the creation of short- and medium-term forecasts for 
the spread of harmful substances, assessment of their impact on the growth of major phytoplankton species in the Azov 
Sea, and determination of strategies for sustainable development management. 
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AHHOTalna 

Beedenue. Matematuueckue MOJeM UM MeTOJbI MOBCEMECTHO HCIOMb3YIOTCA JIA HCCHEOBAHHA TIPHPOAHBIX OOBeK- 
TOB, 34MeHAA Ooee AOporue HaTypHble 9KCTepHMeHTHI. OFHHMU U3 TpyAHOCTeH, BOSHUKAaIOWIMX Ip MOJeHpoBaHHn 
TIPOLECCOB B CJIOKKHBIX CHCTeMaX, ABJIAIOTCA HaJIM4Me BXOAHBIX JaHHbIX MH MOAOop napameTpos Moyen. IIpumenenve 
MCTOJOB YCBOCHHA JaHHbIX HaOIOACHH ABIACTCA OJHUM H3 cIOCOOOB OCHaleHHA MaTeMaTHYeCKHX MOsesel BXOJ- 
HbIMH JaHHbIMH WM 3HaYeHHAMH Mapametpos. Llemb HacTOAMero UcceqOBAHHA COCTONT B MIPOrHO3HpOBaHHH Ha OCHOBe 
MCTOJOB MaTeMaTHYeCKOrO MOJCIMPOBaHHA Pa3sBUTHA COMKHBIX MIPHPOAHBIX CHCTeM B YCIOBHAX 3arpA3HeCHHA Bpel- 
HbIMH BeLecTBaMH. J[If JOCTIDKeHHA eM Ob peWIeHbI cuexyoulue 3aqayH: BbIOpaH MeTOA yCBOeCHHA aHHBIX 
HaOMOCHHH, aKTyaIM3HpoBaHa MaTeMaTH4ecKad MOJeb OMOMOTHYCCKOM KHHETHKH, aHHad MOJelIb CKOMIIJIEKCHpO- 
BaHa C MOJe@JIbIO THAPOAMHAaMHKH, paspadoTaH MporpaMMHbIi Komiiekc. AKTyaJIbHOCTb paOoThI 3akIIOUaAeTCA B IIpH- 
MeHeHHH HOBOTO MOAXONa K pearH3allMH MOJeIM THHAaMHKH (PUTOMIAHKTOHHBIX NOMyIAHH (SBTpopuKayun) A30B- 
CKOFO MOps IPH HaIM4MH 3arpsaSHAIOWIUX IpUMecel, OCHOBaHHOrO Ha IPHMeCHCHHH BapHallMOHHbIX MCTOOB YCBOeCHHA 
JAHHBIX, MOJIYYCHHBIX B XOJ[e KCIeAMIMOHHBIX MCcCeqOBaHHi. 

Mamepuanoi u memoooi. Paciipoctpanenve 3arpa3HAIOWIMX BeELICCTB MOJeuIMpyeTcA Ha OCHOBe TpeXMepHOl MaTema- 
THYCCKOM MO eM, OCHOBaHHOM Ha CHCTeMe ypaBHeHHi KOHBeKUMU — Auddy3nu — peaxuun. Ha Bxoye Moses M0- 
WjaeTCA BEKTOP ABWXKCHHA BOAHOM Cpeybl. CocTaBaloujue BeEKTOpa CKOPOCTH TedeHH B IpHOpexkHOH CucTeMe paccuH- 
TbIBAIOTCA Ha OCHOBE MaTeMAaTHYeCKON MOJeIM THAPOAMHaMUKH, Oa3sHpyrollelica Ha TpexX ypaBHeHHAX JBWKeHHA UM 
ypaBHeHuH Hepa3pbiBHocTH. Pa3spadoTaHHbIi Ha OCHOBe OMMCaHHbIX MOJeseH MporpaMMHbIii KOMIVIeKc MouyyaeT Ha 
BXOle HaTYPHble JaHHble, COOpanHble B XOAe IKCICAMUMOHHBIX HCCIeAOBaHHH, U MO3BOIIACT YTOUHAT MOJCIIb 3arps3- 
HeHHA BOAHOM cpeybI U OMOTEI OarosapsA MPHMeHeHHIO BapHal|MOHHbIX MCTOOB YCBOCHHA JaHHBIX. 

Pesynomamoi uccnedoeanua. WoctpoeH KpaTKOCpOuHbIit MpOrHo3 pacipocTpaHeHHsA 3arpsA3HALOWMX BeLeCCTB Ha BBI- 
xoge u3 Taranporckoro 3amBa. IIpopeqeHHbIii BBIYMCIMTCJIbHbIM IKCICPHMeHT OTpaxkaeT JHHAMHKy pactipocTpaHe- 
HHA 3alpA3HAIOIUMX BELICCTB OT HCTOUHHKOB 3apaxKeHHA Ha BPCEMCHHOM HHTepBalie OT 3 Wo 12 WHe#. 

O6cystcoenue u 3aKii0“enue. PaccMoTpeHHbie B JaHHOM HCCJIeOBaHHH BapHall[MOHHbIe MeTOALI YCBOCHHA JaHHbIxX 
HaOOACHHM WO3BOJIAIOT YTOUHAT HU JOMOMHATL MaTeMaTHYeCKHe MOeIM THHAMUKM (PUTOMIAHKTOHHBIX MOMMyIAWHH 
WM paciipocTpaHeHua 3arpa3HAIOWHx BeulecTB. IIporpamMHoe oGecrieyeHHe, OCHOBaHHOe Ha OMMCaHHBbIX B aHHOl pa- 
OoTe MaTeMaTHYeCKHX MOJeIAX, ACT BO3MOXKHOCTL CTPOUTb KpaTKO- MH CpeHECpOuHbIe MpOrHo3bl paciipoctpaneHusr 
BPeJHBIX IIPHMecel, OLCHMBaTb HX BIMAHMe Ha pa3s3BHTHe OCHOBHBIX BUAOB (PUTOMIAHKTOHHBIX MOMMyIAWHH B A30B- 
CKOM MOpe H OlpeylesIATb CTpaTerHu yupaBsIeHHA YCTOMYMBIM pa3sBHTHeM. 


Ksrouesple cJ10Ba: MO/JICJIb 3BTpoduKallMn, MOJCJIb THAPOTWHaMUKH, BapHalMOHHble MeTOJbI, ONACHbIe ABJICHHA, 
YCBOCHHE JaHHbIX HaOHOeHHii 


Baarogapuoctu. ABtopsi OnaroyqapaT peqakIMOHHYIO KOMaHJy %KYpHasia U AHOHMMHBIX PeCLeH3CHTOB 3a BbICKa3aHHbIe 
3aMe4aHHA, KOTOPbIe MO3BOJIMIIM MOBbICHTb KaYeCTBO CTATbH, a TaKKe BbIPaxKalOT IIPH3HATCJIBHOCTh PYKOBOJUTEIIO 
lpoekta, WieHy-koppecnouyzeuty PAH Anexcanypy Usanosuyy Cyxuuosy. 


@unancuposanne. Vccneqopanne BbINOHeHO 3a cueT rpaHta PoccuiicKoro HayyHoro dona Ne 22—11-00295, 
https://rscf.ru/project/22-11-00295/ 


Aisa uuTMposanua. besoza 1O.B., Huxutana A.B. [pumenenve MeTox0B ycBoeHHa J@HHBIX HaOsFOJeHHM Jy MOJeyIMpoBaHHa 
pactipocTpaHeHHa 3arps3HAIOWIHX BeLeCTB B BOJOeMe Hi YIIpaBIIeHHA YCTOMYMBbIM pa3sBHTHeM. bexonacHocmb mexHOZeHHbIX Uu 
npupooubix cucmem. 2024;8(3):39-48. https://doi.org/10.23947/2541-9129-2024-8-3-39-48 


Introduction. Mathematical models and methods have been successfully used for decades in conducting research in 
various fields of science and engineering. These models provide fast, convenient, and relatively inexpensive tools for 
studying and predicting processes in complex natural systems, compared to expeditions and field experiments. They are 
used to solve a wide range of scientific and practical problems, such as predictive modeling of siltation of shipping 
lanes, which is important for safe navigation, as well as predicting the consequences of emergencies and man-made 
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disasters. An example is the severe storm on November 11, 2007 in the Azov-Black Sea basin, as a result of which 
more than 20 ships were wrecked, and the Kerch Strait area became the site of an environmental disaster. Several tons 
of fuel oil and sulfur got into the water, as a result of which the coastline and sediment layer were contaminated with 
petroleum compounds. Its consequences were observed for several more years. Another example of adverse processes is 
the transport of bottom materials from the mouth of the Don River to the Taganrog Bay, which affects the habitats of 
aquatic organisms, promotes intensive eutrophication and leads to the reproduction of Chironomidae Newman. 

According to Decree of the Government of the Russian Federation No. 2451 dated December 31, 2020', researchers, 
decision makers, and representatives of water protection services have only four hours from the moment of detection or 
from the moment of receipt of information about a spill in a reservoir to calculate changes in the concentration of oil 
and petroleum products and their localization. According to the above and other regulatory documents adopted by the 
Government of the Russian Federation, responsible persons and structures must make a decision and take actions to 
eliminate a dangerous environmental situation of a natural and man-made nature within a few hours or days. Based on 
this, the time for making forecasts and scenarios for the development of an emergency situation is limited. This 
requirement determines the relevance of developing a set of mathematical models of hydrodynamics and hydrobiology 
that would take into account the features of coastal systems (the effect of winds on the structure of currents, the Coriolis 
force, the complex geometry of the computational domain, turbulent exchange, evaporation, wind-driven effects, river 
flows, etc.) and make it possible to obtain forecasts in a minimum time. 

When forecasting the development of natural systems and describing a real physical phenomenon through 
mathematical modeling, it is not enough to simply construct a function of the process state. The use of conjugate 
problems for mathematical models and algorithms based on variational principles can increase the accuracy of 
solutions. These principles allow us to establish a connection between the model and field data [1]. This methodology is 
effective for solving applied problems, including computationally intensive tasks, such as variational and optimization 
problems in mathematical and nuclear physics [2]. G.I. Marchuk and his colleagues used conjugate equations [3], which 
made it possible to improve the efficiency of solving problems of aero- and hydrophysics for the atmosphere and deep- 
water reservoirs. The theory of building conjugate operators for linear and nonlinear models has also been improved [4]. 

A variational approach to solving combined direct and conjugate problems using assimilation methods has improved 
the relationship between mathematical models and field data. Methods of data assimilation, which have been 
developing since the 1960s, are based on the construction of inverse and optimization problems using two approaches: 
the classical Lagrange variational principle using conjugate problems [4] and optimization methods such as weighted 
least squares [5]. 

The assimilation of observational data is a tool that can significantly improve the accuracy of predictive modeling of 
natural processes. It has long been successfully used by the scientific community [6]. Here, the urgent task is to develop 
new methods that would significantly reduce the calculation time. 

Materials and Methods. When building models for predicting natural phenomena and processes, one of the main 
problems was the correspondence of the solution obtained using a mathematical model to the real process occurring in 
the natural system and reducing the percentage of uncertainties. 

When constructing mathematical models of hydrodynamic and hydrobiological processes, information about the 
initial conditions and parameters of the model was required, which could be obtained using observational data. Thus, 
when constructing predictive scenarios of natural or man-made emergencies, it was very important to assess the 
adequacy of the mathematical model itself. The next stage of modeling included checking the correctness and stability 
of the task. The study of a mathematical model at a continuous level involved the study of the influence of input data on 
the solution of a model problem. The perturbation of the right-hand sides of the equation used or the system of partial 
differential equations in the Cauchy problem under consideration made it possible, with known operators, to investigate 
the properties of the constructed mathematical model. The study of stationary and special points of a continuous 
function or several functions — solutions to the task, for example, the concentration of one or more pollutants in the 
aquatic environment, allowed us to develop scenarios, from pessimistic to optimistic ones, in order to develop measures 
for effective management of a complex aquatic ecosystem. 

For the first time, the method of polynomial interpolation of data from a constantly replenished database of 
environmental measurements in the two-dimensional case was used to analyze field data. Areas of influence have been 
identified for observations. When implementing the algorithm for calculating the current background, the data from the 
previously received forecast was used as input information. 


' On Approval of the Rules for the Organization of Measures to Prevent and Eliminate Oil and Petroleum Product Spills on the Territory of the Russian 
Federation, with the Exception of the Internal Sea Waters of the Russian Federation and the Territorial Sea of the Russian Federation. Decree of the 
Government of the Russian Federation No. 2451 dated 31.12.2020. URL: https://docs.cntd.ru/document/5733 19208 (accessed: 21.05.2024). (In Russ.) 
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The OI (Optimal Interpolation) method has brought a new level of sophistication to the methodology for solving 
problems of data assimilation. This method was based on statistical interpolation techniques. 

The next stage in the development of the considered methods was associated with the development and 
implementation of variational methods, including the theory of optimal control. These methods were based on 
minimizing the functionality built in a special way, with the help of which a connection was established between 
solutions and observations (field measurements, expedition data, GIS geographic information system databases). 

This theory and methods were widely used in the implementation of the tasks of meteorology [5] and dynamic 
oceanography [7]. In the process of minimizing the constructed functional, it was necessary to calculate its gradient, for 
which conjugate equations were successfully applied, which was described in [8, 9]. 

Model of the Dynamics of Phytoplankton Populations. A mathematical model of the dynamics of phytoplankton 
populations describes the process of active growth of microalgae in the presence of a sufficient number of biogenic 
elements. If the development of phytoplankton populations becomes too intense, it is referred to as the process of 
eutrophication. The causes of eutrophication can be both natural (climatic changes) and anthropogenic in nature 
(the entry of a significant amount of nutrients into the reservoir from river drains). System C; contains the concentration 
values of the i-th substance [10, 11]: 


(uC O(vC O((w+ we, JC, 
OG @ dy U an (( ct) t) (ut) Ree: Efe blew, (1) 
Ot Ox oy Oz Ox ox oy oy Oz OZ 
where u= {u, v, w} — velocity vector of the medium (water flow); wc, — gravitational deposition of the k-th 


component, if it is in a suspended state; u,, ve — horizontal and vertical components of the coefficient of turbulent 
exchange for the k-th componen; w; — chemical-biological source (drain) or term describing aggregation (sticking- 


splitting), if the corresponding component is a suspension, index k indicates the type of substance, k =1,15: 

1 — hydrogen sulfide (2S); 

2 — element sulfur S; 

3 — thiosulfates (and sulfites); 

4 — sulfates (SOu); 

5 — total organic nitrogen (JN); 

6 — ammonium (NH4 — ammonium nitrogen); 

7 — nitrates (NO3); 

8 — nitrites (NO2); 

9 — phytoplankton; 

10 — zooplankton; 

11 — silicates (SiO; — metasilicate; SiO, — orthosilicate); 

12 — dissolved oxygen (O2); 

13 — ferrum (Fe?*); 

14 — phosphates (POs); 

15 — silicic acid (H2SiO4, H2SiO3 — orthosilicic and metasilicic acids, respectively); 

16 — microplastic. 

System (1) contains equations that can be attributed to the convection — diffusion — reaction type. As a 
computational domain, let us consider enclosed basin G. Undisturbed surface of basin 2, is bounded from above 
by G, X»= (x, vy) — bottom surface from below, co — cylindrical surface limits G from the side. We will use the 


following notation: 2 =X») Uo UL, — piecewise smooth boundary of domain G, time interval 0<t<To. We assume 
that n and un — vector of external normal and is the normal component u to surface =. 


We consider that initial conditions of system (1) are as follows: CiH-0 = Cio(x, y, z), k = 115. 


We combine (1) with the following combined boundary conditions: 


for o: if un<0, then C,=0; if un>0, then <9; for Xo: Ft = 9(C,); at the bottom »;;: Ft --2,Cy, 
IZ IZ 


k =1,16 , we define g as the coefficient of absorption of the k-th impurity by bottom sediments. 

When there is no wind, especially in summer, almost anaerobic conditions may occur in the bottom layers of 
shallow reservoirs, such as, for example, the Azov Sea, the Taganrog Bay, and the Gelendzhik Bay. The reduction of 
water-saturated surface sludge entails the release of iron, phosphates, sulfates, manganese, ammonium and silicates, as 
well as organic compounds into the solution. Complexed models of type (1) (developed by the team of authors the 
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model of hydrodynamics [12]) were used to study the mechanisms of manganese oxidation and reduction, NH4 
assimilation, nitrification, nitrate reduction (denitrification), ammonification, H2S oxidation, sulfate reduction, etc. 
Experiments with model (1) made it possible to study biogenic and oxygen regimes of the coastal system, to analyze the 
mechanism of formation due to anthropogenic eutrophication of marine phenomena of fish and other aquatic organisms. 
Variational Approach for a Three-Dimensional Spatial Mathematical Model of Water Eutrophication. Let us 
write down mathematical model (1) for the computational domain (the Azov Sea) in the form of an operator equation: 


L(CY)= DS +5(C¥)-y-1=0, (2) 


we define C as a vector function of the state of the studied aquatic ecosystem C = {c ,(X,t),k = 1,16} , C=C(x,t)eO(U,), 
Wf, =Gx (0,7, 0) : (x,t) eff,; J (S,Y) — differential nonlinear spatial operator; D — diagonal matrix; 
w= {y k (x,1),k = 1,16} — vector the components of which are source functions; r = {ry (x,1),k = 1,16} — vector the 


components of which contain functions of uncertainties and errors of mathematical model (1) with initial and boundary 
conditions. Dependencies (observation models), coefficients and parameters wc4, Lu, U, V, W, Vs, Input data of initial and boundary 


conditions for model (2), k = 1,16 ; internal parameters of the operators are included in Y € R (i, ). 
Let ¢ = 0, then the initial conditions for (2) take the form: 
C=Ci+h, Y=YV) +6 (3) 


where C? and Y° — a priori estimates of the vector function of the state and the vector of parameters, respectively; 


the uncertainty functions are denoted using € and ¢. 
Let us consider the integral identity: 


I(C,Y,C*)= | E(cy).c )acar =0, (4) 
I, 
here C* — functions conjugate to C(C* € O*(II,)). (4) is a variational formulation of model problem (2), (3), or an 
energy-type functional. Let us rewrite (4) in the following form: 
16 
i(cy.c’)=>" (4c.c’), - [ (ys +n Jciacar 0. (5) 
a Ll, 

The operators of turbulent exchange and transfer are included in terms (AC, C’). 

Hydrodynamics model [12] will be considered as a model of the process. In hydrobiology models, parameterization 
of parameters and the resulting functional dependencies, for example, for describing production and destruction 
processes or the growth of phyto- or zooplankton, will be considered submodels, or observation models. Let us define 
the relationship between measurements and state functions: 

@n =| W(C)] +n(%0), (6) 
where [W(C)],, — vector of submodels (observation models); n(x, 4) — vector of errors and uncertainties; @,, — values 
that we are monitoring. 


m 


Let us define @», on Li)" € L,. In (6), the operation of transferring information from I], to L;” is indicated by 


square brackets. 
Let us expand the modeling system with data from field measurements (we consider them close to accurate), while 
the “quality” functionality will have a look: 


04(€)={(.-[W(C)], J Mxo(0n-[W(C)],)} = ("enn o 
yy 
we define yo as a weight function for determining the configuration of the observation carrier L]”” in L, and the 
integrals over domain L],, representing a measure for (7) in the form C; = Myo(x, 4), where M — weight matrix. 
Let us consider functionals representing generalized characteristics of hydrobiogeocenosis behavior: 
Oy (C)= | Fe(C) te (xt) dGdt = (Fearty)s x <0 (H,). k= ER. 
UW, 
F,(C) — limited and differentiable relative to C € O(I],) functions that we will evaluate. 
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Let us define a functional to minimize uncertainties: 


h 
5} (C)=0} ()+4(n'enn),, (Grr), +((C°-C2)'s(C*-Ch)) + 
: (8) 
(ve -¥2)'C,(¥° -¥!)} sae) /24+I"(GY,C*), £21. 
Re (uy 
We assume that C; — weight matrices, i = 1,4. Let us consider the system: 
a} 
ane = DA,C+J*(C,Y)-w-r=0; (9) 
xh 
ee =(DA,)'C, + A™(C,Y)C; +d, =0; 
oC 
* 0 F, Tt 
Ci (x)|_,, = 0 de = (1 (C)+0.5(n'Cin)): 


C° =C)+Cz'C; (0), t= 0;r(x,t) =Cy'C; (x0); 


0 
a -l ‘ os h Lan 
Y=¥, tC l=! (C,Y,C;); 


t 0 ' 
A(C,Y)C' = a" (C+aC.Y)}] .k21 


a=0 
A\(C’, Y), A; we define as operators of the conjugate problem and derivatives or their discrete approximations in 
time; I’, — functions of the sensitivity of models to changes in parameters; C’= 5C; a — set number. 
Let us consider algorithms for assimilation of data from successive observations coming from various observational 
tools into a real-time modeling system. To do this, we will use the methods of splitting and decomposition: 


N,-1 N,-l p 
P= Ot Hh, = 6" x[t tn O(CCLY.e)= >) > oh, (10) 
n=1 n=l] /=1 
here @!, — part of functional (9) for ([t1,t,] on the l-th stage of _ splitting, 


n=1,N,,p we define as the total number of splitting steps. We will carry out the discretization on the basis of additive- 


averaged splitting schemes. We will use an algorithm to find a solution in domain I! with regular uniform time grid 


o' = {tq.n =0,N i} . For research, we will use a structure with phase spaces: 


Pp 
{cr, C7", r7.1=1p} =| Joh (ut) co* (ui). 
i=l 
Method of Solving the Problem. For solution (9), we consider that n =1,N, . Let us describe the algorithm of the 


method used step by step. 
1. Switch to the subgrid decomposition structure, while t = ¢,,_,: 


Pp 


fen" <o"(u)}. Ufer' <ot(u)}, er <en 


i=l 
2. In the subgrid structure, obtain solutions to direct and conjugate problems: 
A'C} -w} -r! =0, 1=1p, p21. 


n-l 


OD, (C 
APC" = [Pel sure (0, -[w(c)], ) a 
I 


* ale sé 
Cr =0, 7 =(CZ) Cy", t, <tst,. 
Functions y7' include C7"'. rj’ takes into account all the uncertainties in the step [ty1, tu]. 


3. Return to main structure QO" (11") att=th 
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P P 
n n ui n 1 n 

Utci <a?(ur)} = {e" eo" (ur)}, c= dic. 

[=1 [=1 

We believe that the last stage of splitting can be implemented using formula: 


A,” =Wp-k, =0; (11) 
here A,,C” — approximation operator part of the model at the p-th stage; w’, — source functions; r;, — function of 
uncertainties of model (2) and the splits introduced into the discrete model at a step. We believe that matrices of 
weights C), and C2, are known. 

Let us consider the following task: 
Nyy C™" = O1y,Cy, ("1 -C""); (12) 


pn 
n_ -l *n 

F, =(C;, (ay, |C™ (13) 
For the solution, we will use the following method, the algorithm of which we will describe step by step: 
1. Calculate C*” using (12). 
2. Find r; , using (13). 
3. Find C’, using formula (11). 
The problem that arises during the implementation of the algorithm of the system of linear equations method will be 


solved by the sweep method. 
Let us consider the transformed algorithm: 


AS,C™ = O4,Cin(@" -C"), (14) 
A pC" —" - (Cz, / 09, )C* =0. (15) 


We will solve the system by the sweep method. For the modifications considered, the stability of the splitting 
schemes ensures the stability of the data assimilation schemes used. If in (8) C* = const, then it gives a balance ratio of 
the first order. If C” = C, then we obtain the equation of the energy balance of the analyzed system. 

Let us write out the quality functionality for the new subsequent modification: 


®,, (C) = 0,5] a, (niin) + oo (15Maut) |. (16) 
r, = Rags Wns Mn = 9), =O), OQ) +A) = L 1), > 0. (17) 
We find the minimum of the functional relative to function C”: 
* n n a n n n 
A'nCin (ApnC” — 8 )+—" Cn (C"-@")=0. (18) 
2n 


We obtain a uniquely solvable system with a five-diagonal matrix. Let us solve the resulting system of linear 
equations using the sweep method. 

A new class of real-time assimilation methods includes a scheme of additive sequential assimilation. Due to the 
large amount of computational work, the built modifications are focused on supercomputing systems, including cluster 
systems and graphics accelerators. 


Measurements @,, are used in the form of maps and digital images. This representation gives a significant density 


of data in domain I], measurements are information fields. The planning of observations is based on the values of the 
uncertainty function. If these values are high, additional field measurements or observations are scheduled. 

Results. To solve the problem of modeling the eutrophication of the Azov Sea waters (1), a set of parallel programs 
has been developed, including: 

—a module of hydrodynamic processes that calculates the flow field of a water stream based on a mathematical 
model for a shallow reservoir [12]; 

—a module for the spread of pollutants in the aquatic environment and changes in the concentration of basic aquatic 
organisms (1), which allows us to assess the effect of pollutants on the biological productivity of the water area; 

—a map of depths of the Azov Sea for the construction of computational grids for the numerical implementation of 
the developed algorithms; 

—a database of expedition data, which allows refining the model of pollution of the aquatic environment and the 
spread of biota through the use of methods of data assimilation described above. 
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The numerical experiment was conducted using the developed software. The velocity vector of the water flow 
was calculated using a hydrodynamic model with an easterly wind at a speed of 5 m/s and fed into the input 
for calculating the movement of pollutants containing microplastic particles based on convection-diffusion equations. 
Figure 1 a shows the results of the numerical experiment to calculate water flow current fields under specified 
meteorological conditions. Vortex structures of currents were observed in the region of spits, in the northeastern part 
of the sea in Taganrog Bay. The color gradient in Figure | a represents the distribution of water flow velocity, with 
the maximum value at 4.822 m/s. Figures | b, c, and d show the results of calculations for a hypothetical scenario of 
the spread of microplastic contamination in the aquatic environment, considering the input of microplastics into the 
Azov Sea through the drains of the Don and Kuban rivers and the presence of a source of hazardous substances 
including microplastics at the outlet of Taganrog Bay. The figures also include a short-term forecast for the spread of 
contaminants 3, 6, and 12 days after the initial input. The initial concentration of pollutants was 5 mg/L, and after 
3 days, the maximum concentration at the Taganrog Bay outlet was 1.363 mg/L. After 6 days, it decreased 


to 0.83 mg/L and after 12 days, it dropped to 0.336 mg/L. 
- 1.363 
mg/l 
1.022 
4.822 
2.893 
1.929 0.341 
0.964 
0.000 0.000 
b) 
mg/l 
0.750 
0.500 
0.250 
0.000 
d) 


a) 
Fig. 1. Currents in the Azov Sea and the spread of pollutants, time interval: 


°) 


a — initial concentration; b — 3 days; c — 6 days; d— 12 days 


Discussion and Conclusion. The computational experiment shows that, despite the easterly wind, vortex structures 
in the currents captured pollutants and transported them to the Taganrog Bay. Stable vortices have the potential to 
capture and retain microplastic particles that enter the sea from river drains. They also contribute to the accumulation of 
pollutants in the lower layers due to biofouling by microplastic particles and their sinking. 

As mentioned above, when creating mathematical models to predict natural phenomena and processes, one of the 
main challenges is to validate their accuracy by analyzing the outcomes obtained from them to ensure they match the 
behavior of the natural system under study. When creating mathematical models of hydrodynamic and hydrobiological 
processes, information on initial conditions and parameters (input data) is needed, which can be acquired through 
observations. Therefore, when creating predictive scenarios, it is essential not only to assess the quality of the 
developed mathematical model but also to incorporate observational data and investigate the sensitivity of these models 
to variations in input data. 

The paper presents an approach to implementing a model of phytoplankton population dynamics (eutrophication) in 
the Azov Sea using variational methods for assimilating data obtained during field studies. The software package 
developed uses materials from the field studies, constantly updated environmental databases, and GIS, and allows 
refining the model of aquatic pollution and the spread of aquatic organisms using variational data assimilation 
techniques. The developed software makes it possible to forecast the spread of pollutants in the coastal system, 
including biogenic substances that act as a nutrient medium for the growth of dangerous microalgae, leading to 
eutrophication. This forecast enables the development of strategies for sustainable management of the natural system 
and its protection. 
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Sanenennoilit eknad aesmopoe: 

¥O.B. Bestopa: pa3padoTka MmporpaMMHOTO MO/JyIA, MpoBeyeHwe IporHocTM4ecKoro pacueTa JBYXKeHHA 
3alpA3HAFOIIMX BelecTB, OPOpMsIeHve Hay4HOW CTaTBH. 

A.B. Hukutwna: onvMcanve TeopeTHueckol 4aCTH HCCIeqOBaHHA pacliIpoOocTpaHeHHA 3arpA3HAIOWIMX BeIeCTB C 
IIPHMeHeHHeM BapHaliMOHHBIX MeCTOJOB YCBOCHHA JaHHBIX HaOIOTeHuH. 
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